*     exponential integral
      DOUBLE PRECISION FUNCTION RFE1(DU)
      IMPLICIT NONE
      DOUBLE PRECISION DU

      DOUBLE PRECISION DA1, DA2, DA3, DA4, DA5, DA6,
     &                 DAF1, DAF2, DBF1, DBF2
      PARAMETER (DA1=-.57721566, DA2=.99999193 , DA3=-.24991055,
     &           DA4=.05519968 , DA5=-.00976004, DA6=.00107857 )
      PARAMETER (DAF1=.250621  , DAF2=2.334733)
      PARAMETER (DBF1=1.681534 , DBF2=3.330657)

      DOUBLE PRECISION DE1

      IF (DU .LE. 1.) THEN
         DE1=((((DA6*DU+DA5)*DU+DA4)*DU+DA3)*DU+DA2)*DU+DA1
         DE1=DE1-LOG(DU)
      ELSE IF (DU.LE. 111.) THEN
*        prevent overflow in DEXP (111. large enough not to lose accuracy
*        after the conversion to single precesion: DEXP(111)=1.609487e+48)
         DE1=((DU+DAF2)*DU+DAF1)/((DU+DBF2)*DU+DBF1)
         DE1=DE1/(EXP(DU)*DU)
      ELSE
         DE1=0.
      END IF
      RFE1=DE1

      RETURN
      END
